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Abstract: The positive weight next-to-leading-order (NLO) matching scheme (POWHEG) 
is applied to Drell-Yan vector boson production in the Herwig++ Monte Carlo event gen- 
erator. This approach consistently combines the NLO calculation and parton shower sim- 
ulation, without the production of negative weight events. The simulation includes a full 
implementation of the truncated shower required to correctly describe soft emissions in an 
angular-ordered parton shower, for the first time. The results are compared with Tevatron 
and Z production data and predictions for the transverse momentum spectrum of 
gauge bosons at the LHC presented. 
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1. Introduction 

Monte Carlo simulations traditionally combine leading-order 2^2 matrix elements with 
parton showers which provide resummation of soft and collinear radiation. This provides a 
fully exclusive description of observables which allows evolution to the hadronization scale 
where models of the non-perturbative regime can be incorporated. This makes Monte Carlo 
simulations an essential tool in experimental analysis allowing a fully exclusive description 
of the final state to be compared directly with experimental results. 

Parton shower simulations include the leading-logarithmically, and an important sub- 
set of the next-to-leading-logarithmically, enhanced contributions and therefore underesti- 
mate the radiation of high transverse momentum (pt) partons. The Monte Carlo description 
can be improved by matching to higher-order matrix elements which correctly describe the 
production of high px particles. A number of approaches have been developed to correct 
the emission of the hardest parton in an event. In the PYTHIA event generator [1], correc- 
tions were included for e + e~ annihilation [2], deep inelastic scattering [3], heavy particle 
decays [4] and vector boson production in hadron collisions [5]. In the HERWIG event gen- 
erator [6,7] corrections were included for e + e~ annihilation [8], deep inelastic scattering [9], 
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top quark decays [10] and vector boson production [11] in hadron-hadron collisions follow- 
ing the general method described in [12,1s]. 1 These corrections had to be calculated for 
each individual process and only corrected the first or hardest 2 emission, in addition the 
method can only be applied to relatively simple cases and the leading-order normalisation 
of the cross section is retained. 

In recent years there have been a number of additional developments which aim to 
improve on these results by either providing a description of the hardest emission together 
with a next-to-leading order (NLO) cross section [19-31], 3 or the emission of more than 
one hard parton at leading order [32-37]. 4 These matching prescriptions are complicated 
because the regions of phase space filled by the higher-order matrix elements and the parton 
shower must be smoothly separated in order to avoid problems such as double-counting, 
where the shower and matrix elements radiate in the same region. In general the major 
complication is gaining an analytic understanding of the result of the parton shower either, 
to subtract it from the real emission matrix element, as in the MC@NLO approach [19-24], 
or to reweight the real emission matrix elements so they can be merged with the parton 
shower in multi-parton matching. 

The first successful scheme for matching at NLO was the MC@NLO approach [19-24] 
which has been implemented with the HERWIG event generator for many processes. The 
method has two draw backs; first, it involves the addition of correction terms that are not 
positive definite and therefore can result in events with a negative weight and second, the 
implementation of the method is heavily dependent on the details of the parton shower 
algorithm used by the event generator. 

In Ref. [26] a novel method, referred to as POWHEG (POsitive Weight Hardest Emis- 
sion Generator), was introduced to achieve the same aims as M CO NLO while creating only 
positive weight events and being independent of the event generator with which it is im- 
plemented. The POWHEG method has been applied to Z pair hadroproduction [27], heavy 
flavour hadroproduction [30] and e + e~ annihilation to hadrons [31]. A general outline of 
the ingredients required for POWHEG with two popular NLO subtraction schemes is given 
in Ref. [29]. 

The POWHEG shower algorithm involves generating the hardest emission in pr sepa- 
rately using a Sudakov form factor containing the full matrix element for the emission of an 
extra parton and adding to this vetoed showers, which produces radiation at lower scales in 
the shower evolution variable, and a truncated shower, which generates radiation at higher 
scales in the shower evolution variable, than the scale of the highest pt emission. While the 
POWHEG scheme is independent of the parton shower algorithm, it does require the parton 
shower to be able to produce vetoed and truncated showers. The ability to perform vetoed 

1 Herwig-|— I- includes matrix element corrections for e + e~ annihilation [14], top quark decays [15], vector 
boson production [16] in hadron-hadron collisions and Higgs production in gluon fusion [17] using the 
approach of Refs. [12, 13] and the improved angular-ordered parton shower of Ref. [18]. 

2 In PYTHIA the first emission was corrected whereas in HERWIG any emission which could be the hardest 
was corrected. 

3 There have been other theoretical ideas but only the MCONLO and POWHEG methods have led to 
practical programs whose results can be compared with experimental data. 
4 A recent comparison of these approaches can be found in [38]. 
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showers is present in most modern Monte Carlo event generators, however some changes are 
required to enable them to generate truncated showers. Although the POWHEG approach 
is formally correct to the same accuracy as the MC@NLO technique the two methods differ 
in their treatment of sub-leading terms. 

In this work the POWHEG approach is applied Drell-Yan vector boson production with 
the Herwig++ [17,39] event generator. A full truncated shower is implemented for the first 
time. 5 

The paper is organised as follows. In Sect. [2] we outline the main features of the 
POWHEG method. In Sect. ^ we collect the essential formulae relating to the NLO cross 
section, for implementation in the program. In Section [4] we give details of the event 
generation process for the hard configurations, this is followed by a description of how 
these configurations are subsequently reproduced by the Herwig++ angular- ordered parton 
shower, thereby accounting for colour coherence effects associated with soft wide angle 
parton emissions. In Sect. [5] we present the results of our implementation, comparing it to 
Tevatron data, and in Sect. |^ we give our conclusions. 

2. The POWHEG method 

In the POWHEG approach [26] the NLO differential cross section for a given iV-body process 
can be written as 

da = B d$B 
where B ($b) is defined as 

B ($ B ) = B + V (* B ) + / R (* B , *r) ~J2 Ci d$ «' ( 2 - 2 ) 

J i 

B ($b) is the leading-order contribution, dependent on the iV-body phase space variables 
the Born variables. The regularized virtual term V (3\b) is a finite contribution arising 
from the combination of unresolvable, real emission and virtual loop contributions. The 
remaining terms in square brackets are due to iV+1-body real emission processes which 
depend on both the Born variables and additional radiative variables, <&r, parametrizing 
the emission of the extra parton. The real emission term, R ($_b, $r), is given by a sum of 
parton flux factors multiplied by real emission matrix elements for each channel contributing 
to the NLO cross section. Finally, each term C, (&b,&r) corresponds to a combination of 
real counterterms / counter- event weights, regulating the singularities in R ($>b, &r)- The 
modified Sudakov form factor is defined as 

& R ipr) = exp 



6 Truncated shower effects were neglected in Refs. [27,30], while an approximate treatment was used in 
Ref. [31] where at most one truncated emission was generated. 
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where kj< &r) is equal to the transverse momentum of the extra parton in the collinear 
and soft limits. 

In the framework of a conventional parton shower Monte Carlo program the Sudakov 



form factor A^, (jyp) has the same form as that in Eq. 2J3 but with R(&b,®r) replacing 
R(<&Bi&r), where the former is typically the sum of the real emission matrix elements 
approximated by their soft /collinear limits. This object has an interpretation as the prob- 
ability that given some initial configuration of partons resolved at some characteristic scale 
p T , in evolving to some other resolution scale, pt, no further partons are resolved. The 
parton shower approximation to the NLO differential cross section is analogously given by 



da = B d$ B 



A « (0) + bW) Ar {hT d * R 



(2.4) 



The first term in Eq. gives a differential distribution proportional to the leading-order 
differential cross section, which is retained in the absence of any parton showering, i.e. with 
probability A^ (0) (non-radiative events), while the second term represents the probability 
of evolving from the starting scale to scale pt at which point an emission is generated 
according to R ($b, &r) ■ The term in square brackets in Eqs. 2.1 and 2.4 is equal to 

PT 

unity when integrated over the radiative phase space. 

In the POWHEG framework the parton shower is promoted to NLO accuracy by de- 



manding that the non-radiative events are distributed according to the first term in Eq.|2.1 
and that the hardest (highest px) emission is distributed according to the second term. 
Whereas in the conventional simulation an iV-body configuration is generated according 
to B (<3?_b) and then showered using the Sudakov form factor A^, the POWHEG formalism 
requires that the iV-body configuration is generated according to B (&b) and then showered 
with the modified Sudakov form factor A^. As the B ($b) term is simply the NLO differ- 
ential cross section integrated over the radiative variables, it is naturally positive, which is 
reflected in the absence of events with negative weights, which are an unpleasant feature of 
other next-to-leading order matching schemes. 

Since any further emissions constitute higher-order terms in the differential cross sec- 
tion, next-to-next-to-leading order and beyond, we may generate higher multiplicities in the 
usual way, by showering the radiative events using the standard parton shower algorithm. 
These showers must not generate emissions with pr greater than that of the emitted parton 
in the event generated according to Eq. 2.1, 6 they must be vetoed showers. For a parton 
shower which evolves in pt the POWHEG implementation is a trivial task; one simply ini- 
tiates a parton shower from each external leg of the radiative POWHEG event using its px 
as the initial-evolution scale. 

Angular-ordered parton showers account for the phenomenon of QCD coherence where 
wide-angle soft gluon emissions, from near collinear configurations of two or more partons, 
have insufficient transverse resolving power to be sensitive to the constituent emitters. In 
effect the resulting radiation pattern is determined by the colour charge and momentum 



Should the POWHEG event turn out to be a non-radiative event, the pr veto scale is zero so these 
events remain as no-emission events. 
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of the mother of the emitters, rather than the emitters themselves. Ordering branchings 
in the parton shower in terms of the angle between the branching products takes this 
logarithmically enhanced effect into account; wide-angle soft partons are emitted from the 
mother of any subsequent smaller angle splittings by construction. 

As well as circumventing the problem of negative event weights, a second major suc- 
cess of the POWHEG method is in defining how the highest pt emission may be modified to 
include the logarithmically enhanced effects of this soft wide-angle radiation. In Ref. [26] 
it was shown how the angular- ordered parton shower which produces the hardest emis- 
sion, may be decomposed into a truncated shower simulating coherent, soft wide-angle 
emissions, followed by the highest px (hardest) emission, followed again by further vetoed 
parton showers, comprising of lower pr, smaller angle emissions. Performing this decompo- 
sition established the form of the truncated and vetoed showers, thereby describing all of 
the ingredients necessary to shower the radiative events in the POWHEG approach. This 
procedure was proven in [29] to give agreement with the NLO cross section, for all inclusive 
observables, while retaining the logarithmic accuracy of the shower. 

In the POWHEG framework positive weight events distributed with NLO accuracy can 
be showered to resum further logarithmically enhanced corrections by: 



generating an event according to Eq. 2.1 



• directly hadronizing non-radiative events; 

• mapping the radiative variables parametrizing the emission into the evolution scale, 
momentum fraction and azimuthal angle (q^, Zh, 4>h), from which the parton shower 
would reconstruct identical momenta; 

• using the original leading-order configuration from B ($b) evolve the leg emitting the 
extra radiation from the default initial scale down to q~h using the truncated shower; 

• inserting a branching with parameters (q~h, Zh, (f)h) into the shower when the evolution 
scale reaches q\\ 

• generating pt vetoed showers from all external legs. 

This procedure allows the generation of the truncated shower for the first time with only a 
few changes to the normal Herwig++ shower algorithm. 



3. Next-to-Leading Order Cross Section 

Although the NLO cross section for the Drell-Yan process was calculated nearly 30 years 
ago [40,41], we have implemented an independent calculation of it more suited to our present 
goal, including the decay of the vector boson and j/Z interference effects. In this section 
we collect the ingredients that arise in the NLO calculation for q + q — ► I + I, necessary to 
describe our implementation of the POWHEG method. 
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3.1 Kinematics and phase space 

The leading-order process under study is of the type, p® + Pq — > Pi + ■•■ +Pn, in which 
all particles in the iV-body final state are colourless. We denote the incoming hadron 
momenta P@, for hadrons incident in the ±z directions, respectively. The corresponding 
massless parton momenta, with momentum fractions and Xq, are given by p@ = x @ P @ . 
The momenta of the particles produced in the leading-order iV-body process are pi, where 
i ranges from 1 to N. 

We use $s to denote a set of variables defining a point in the iV-body phase space and 
convolution over the incoming momentum fractions, and $b to denote a set of variables 
parametrizing the iV-body phase space in the centre-of-mass frame of the partonic process. 
It will also be convenient to define p as the total momentum of the colour neutral particles, 
p = x®P© + XqPq, and y as the rapidity of p. The partons' momentum fractions are 
therefore 




s 

The phase space for the leading-order process is 




x e = \-e~y. (3.1) 



d$ B = dx e dx e d6 B = -dp 2 dydi B , (3.2) 

s 

where d<&B is the Lorentz invariant phase space for the partonic process, in n = 4 — 2e 
dimensions, 

dl> B = (27r) B [J ^ Sn (Pffi + Pe ~ P) ■ m = ,„ , n -i^ . (3-3) 

and s the hadronic centre-of-mass energy. The set of variables &b = |p 2 , y, defines 
the Born variables. 

The real emission corrections to the leading-order process consist of 2 — > N + 1 
processes, + Pe — > pi + ... +pat + fc, where we denote the momenta of the N final-state 
colourless particles pi and that of the extra colour charged parton by k. The momentum 
fractions of the incoming partons are distinguished from those in the 2 — > N process as Xq, 
and Xq (p@ = x @ P@). For these processes we introduce the Mandelstam variables s, t, u 
and the related radiative variables $i{ = {x, v, (f>}, which parametrize the extra emission: 

s = (P®+Pef = — > (3-4a) 

2 

i = (Pe -kf = P -{*- 1) (1 - v) , (3.4b) 

u = (pe — ^) = — (x — l)v , (3.4c) 

where 4> is the azimuthal angle of k with respect to the +z axis, and p the total momentum 
of the colourless particles, p = x^Pq + XqPq — k. 
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To perform a simultaneous Monte Carlo sampling of the N- and iV+l-body phase 
spaces one has to specify the integration variables. We choose two of these to be the mass 
and rapidity of the system of colourless particles, therefore p 2 = p 2 and y = y, where y is 
defined by analogy to y as the rapidity of p. 7 The momentum fractions of the partons for 
2 — ► N + 1 processes are therefore related to those of the 2 — > N process by 




*e = -W 1 r (1 Tf )t; v (3-5) 



l-(l-x)v ' w y/x]l l-(l-x)(l 

The phase space of the iV+l-body real emission processes can be written in n = 4 — 2e 
dimensions as 

V 2 /4ttY 1 
d** +1 = dd> B d** -fjj-j ^_ j fir — y J" (x, „) , (3.6) 

where here the partonic Born variables specify a configuration in the rest frame of p 
rather than p. The function J (x, v) is given by 

J (x, v) = [S5 (1 - x) + C (x) (5 (v) +S(l-v))+H (x, v)} v (1 - v) (1 - xf , (3.7) 

where 

1 ^2 

(3.8a) 

lnx + 2[ 1 ^-^-) , (3.8b) 



s = 


1 


vr 2 
6 ' 


C(x) = 


1 

e 


1 

(1-x 


H (x, v) = 






(1- 


-x) + 1 



- (!-«)+" V 1 — a? / + 
1 1 



w+ (i - «; 



+ 



(3.8c) 



The labelling S, C,H reflects the fact that the S and C terms are multiplied by (5-functions 
which limit their contributions to configurations with soft (x — > 1) and collinear (v — > 0, 1) 
emissions, while W is not associated with soft or collinear configurations but instead con- 
tributes to hard emissions 8 of the extra parton k. 

The radiative phase space can be parametrized in terms of the radiative variables 

d$ R = — dx dv d(f> , (3.9) 

where, in the partonic centre-of-mass frame, x = 1 — k°/E, with E the energy of either of 
the colliding partons, and v = \ (1 + cos#), with and <f> the polar and azimuthal angles 
of k with respect to the +z axis. 

As the rapidity of p andp are equal, it is always possible to define a boost B = B^ 1 Bt Ml, 
such that lp = p, where Ml is a longitudinal boost to the frame in which y = and Mt is 
a boost in the transverse direction, to the frame in which the transverse momentum of p is 
zero. It follows that an iV+l-body configuration can be assembled by first reconstructing 



7 Henceforth we will always refer to these variables as p and y. 

s Here by hard we simply mean emissions which are neither soft or collinear. 
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the iV-body configuration corresponding to &b, then p and k (from p 2 , y, <£_r), at which 
point the boost B can be calculated and its inverse applied to the Af-body configuration. 
Although Ml is uniquely defined due to p and p having the same rapidity, y, the transverse 
boost By can be modified according to By — ► WM?, with M a rotation, and B will still 
satisfy Mp = p. A convention must be adopted to fix By, we shall return to this point later. 

3.2 Matrix elements 

The squared, spin and colour averaged leading-order matrix element for the q + q — ► I + I 
Drell-Yan process is given by Mqq (p q ,Pq) , where the first (second) argument refers to the 
incoming fermion (antifermion) momentum. The real emission radiative corrections consist 
of three processes: q + q—>l + l + g; g + q — > I + 1 + q; and g + q — > l + l + q. The matrix 
elements squared for these processes are given by 



M 



R 

qq 



qq 



p 2 tu 



(S + t) Mq q (p q ,pg g ) + {S + U) Mqg(p qg ,Pq) 

-\e{t + U) 2 (Mqq{Pg,Pgg) + Mq q (Pgg,Pg)) 



M 



R 

q9 



q9 



p 2 us 



M% = - 



p 2 st 



(t + u) M q q (Pqg,Pq) + (t + s) M qq - (p~ q ,Pqg) 

(tt + S) 2 (M^ (Pg,Pgg) + M^g (p 9 g,Pg)) 
(U + sf Mqq (Pqg,Pq) + {u + t) 2 Mqq (p~ q ,Pqg) 

~\ e {* + ^ ( M qq (Pg'Pgg) + M qq 



(3.10a) 



(3.10b) 



(3.10c) 



where, for a more uniform notation, and to help show how the crossing of the leading-order 
process is manifest, we have denoted the final-state quark momentum in the qg initiated 
process by pq and the final-state antiquark momentum in the gq process by p q . The shifted 
momenta pi, pj g , and the normalization constants N q q, M qg are given by 



Pi = —Pi ■ 

Xi 



Nqq = &TTasCFH 2e , 



Xi — 



2p.pi 



M qg = STrasTp/j 26 / (1 - e) , 



(3.11) 



Pjg =P~Pi: 



where \x is the regularization scale emerging from the use of conventional dimensional reg- 
ularization. The shifted momenta satisfy 



Pi 



0, 



Pi + Pjg = P, 



(3.12) 



i.e. they obey the relations required for them to be considered as describing a kinematic 



configuration for the leading-order process, hence they form valid arguments for M qq -. 

These matrix elements hold independently of the type of exchanged vector boson, 
they also hold for the case that the leading-order process consists of interferences between 
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diagrams with different exchanged vector bosons. This universal behaviour is due to the 
factorisation of the NLO hadron tensor, into kinematic factors multiplying the leading- 
order hadron tensor. Such a factorisation of the matrix element is not necessary for the 
implementation of the POWHEG method but it improves the flexibility and generality of our 
implementation of the Drell-Yan process. This allows us generate iV+l-body configurations, 
according to the full real emission matrix element, given the leading-order configuration 
using a technique of sampling the radiative phase space with a branching algorithm known 
as the Kleiss trick [13,42], which we have extended to n dimensions as needed for a complete 
NLO calculation, as we now describe in detail. 

The O (as) virtual corrections to the q + q — ► l + l process solely consist of the qq vertex 
correction. At NLO, this loop diagram only contributes to the matrix element through its 
interference with the leading-order amplitude, correcting it by, 



M 



C*sCf ( 47Tjti 



2tt 



1 



r(i 



3 
e 



8 + 



(3.13) 



3.3 Differential cross section 

The partonic flux due to parton a in hadron A and parton b in hadron B, at scale fi 2 , with 
momentum fractions and x e respectively, is defined as 



(3.14) 



In C a b (x^, Xq) the functions // (xj,// 2 ) are the parton distribution functions (PDFs) for 
finding a parton i in hadron / with momentum fraction x\ at scale fi 2 . The contribution to 
the differential cross section from the leading-order process q + q — > I + I is therefore 



daf q = B (* B ) d$ B 



(3.15) 



where 



The virtual corrections (Eq. |3.13|) add 



(3.16) 



where 



V (* B ) 



2tt 



with e defined by - 



l _ l 



dal- = V ($ B ) d$ B , 



2/ri +ln^ 



7b + In (47r) and 



3 ^ „ ' 



B 



(3.17) 
(3.18) 



-31n(^U 2;i " 



(3.19) 



In Eqs.^lH and |T7 we have included a flux factor l/2j? 2 for the partonic process. The 
subscript on Vq identifies the bare divergent quantity. 
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The differential cross section for the real emission processes a + 6^/ + / + cisof the 
form 



d<S> B d<S> R , (3.20) 



where for i = qg, q, g, gg, we have j (i) = q, qg, gg, g. Each of the fi (x, v) functions 

<R 

l ab 



is defined as the coefficient of the squared leading-order matrix element in the M^b rea l" 



emission matrix element, Eq. 3.1C 



Returning to the earlier discussion of the definition of the Born variables and the boost 
B, we note that for each i in Eq. |3.20| the boost to the rest frame of p, + pj is the same up 
to an overall rotation (as pi +f)j = p). Furthermore, for i = q we have pi = Pi/xi, so pi and 
Pi have the same orientation, as do pj and pj for index j = q. Since B was defined up to an 
arbitrary rotation M, we could choose to resolve this ambiguity by setting E = R g , B = B g , 
where M q additionally satisfies, 

Bgp g =Bl 1 R g B T M L pg=pQ, (3.21) 

i.e. Eg is such that the value of p q in the p rest frame is equal to j% in the p rest frame, 
therefore 

M q 3 q -(p q ,Pq g )d$B = M q 3 q -(p (3 ,Pe)d$ B . (3.22) 

To have the analogous equivalence for the i, j = qg, q term, one needs to have B = M qg , 
defined such that ^> qg Pq = Pq is satisfied along with original requirement (Bp = p). 

Similar boosts can be constructed for the i = g, gg terms in the matrix element, 
mapping p g to p® and p© respectively. However, as these terms carry a factor e, they only 
contribute to exactly collinear configurations in which the parton k is unresolved, hence 
those boosts are a purely formal consideration, and are not needed in practice. 

The real emission phase space can therefore be sampled using a simple Monte Carlo 
branching algorithm; namely, given a set of Born and radiative variables, the event is re- 



constructed as described in Sect. 3J_, except that, where previously a single boost B _1 was 
used to embed the iV-body configuration in the iV+l-body event, now we use a boost 
B^ 1 selected according to the probability distribution Vi = fi(x,v) / Yli fi{ x i v )- Sam- 
pling the phase space in this way, the generation of the Born variables becomes completely 
independent of the radiative variables: 

daH = R ab>0 (* B , d$ R d<Z> B (3.23a) 

s - i— * d$R 

Rab,o (*b, *r) = ^ Cab (x e ,x e ) fi 0> v ) ^ B (®b) , (3.23b) 

where C a b is the ratio of the general parton flux relative to that of the leading-order process, 

' ( Ie , Ie) -§±|5a§8> (3.24) 



The functions R a b,o (&b,&r), calculated according to Eq. ^.23| , are given by 

Rab.O {$B, *b) = - Lb (x e , x e ) TZ abfi B (* B ) , (3.25) 

Z7T X 
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where C a h = Cf for ab = qq and Tp otherwise. 
For the qq contribution IZqqfl is given by 

Kqqfl = SqqS (1 - x) + [~Pqq + Cqq^j (5 (v) + 8 (1 - «)) + W,,, (3.26) 

where 

S ^(? + 3)(l +1 n(^))4- 31ll ^), 

P M = (i±^) + , (8.27b) 

c - - c + * 2 > (o4; ln (£) + 2 Gt^t) J + 1 " 1 • (3 - 27c) 

** = (T^ + (T^) (d - -) 2 (1 - 2» (1 - »)) + 2.) . (3.27d) 
For the contribution 7£,j ff) o is given by 

Kqgfi = (-| + C M ) 5 («) + H qg , (3.28) 



where 

P OT = x 2 + (1 - x) 2 , (3.29a) 
C g9 = (x 2 + (1 - x) 2 ) ^ln ^-^-^ + 2 ln (1 - x)^j + 2x (1 - x) , (3.29b) 

ft M = — (2x (1 - x) v + (1 - x) 2 v 2 + x 2 + (1 - x) 2 ) . (3.29c) 

The function H g qfl is equal to 1Z qg fl under the replacement v <-> 1 — v. 

The soft term in the real correction R qq fl, proportional to 5(1 — x), is now combined 
with the virtual correction Vo, to give a soft-virtual contribution V, which is finite for e — > 0. 
The remaining divergences are initial-state collinear divergences proportional to 5 (v) and/or 
5(1 — v). Working in the MS scheme these are absorbed into the definition of the PDFs. 
This renormalization and cancellation of divergences amounts to the replacements, Vo — ► V 
and R abfl -> R ab , where 

cxsCf 



and 



V{$B) = ^r JL V qq -B(<S> B ) (3.30) 

Rab ($B, $/?) = " ^a6 4ft (*©, Zq) 5 ($ B ) , (3.31a) 
Z7T X 

^ ($b, = C OT - (^ («) + <5 (1 - «)) + (3.31b) 

ft M ($b, = (v) + (3.31c) 

(^^' = C 9? 5 (!-«)+ ^ot- (3.31d) 
Adding these contributions we obtain the NLO differential cross section: 

drr = B ($ B ) &§ b + V ($ B ) d$ B + i? ($ B , d$ B d$fl , (3.32) 
where i? denotes the sum of the R ab terms. 
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4. Implementation in HerwigH — \- 



In the first two parts of this section we describe how distributions of NLO accurate non- 



radiative and single emission events are generated. In Sect. [4^ we describe the simulation 
of further, lower pr, emissions, from the radiative events, using the truncated and vetoed 
shower algorithms. 

4.1 Generation of the leading-order configuration 

As noted in Sect. ^, sampling the B(&b) function (Eq. ^.2| ), which is the next-to-leading 
order differential cross section integrated over the radiative variables, 



B($ B ) =B(* b ) 



1 + V<? 9 + / d ® R \ nah ^ ah x e) 

ab 



(4.1) 



provides Born variables distributed according to the exact NLO differential cross section. 
The way in which the leading-order process is factorised inside the real emission terms R a b 
allows the B ($b) distribution to be generated as a straightforward reweighting 9 of the 
leading-order cross section. 

For convenience the radiative phase space &<&r is reparametrized by variables on the 
interval [0, 1] such that the radiative phase space volume is unity, a three-dimensional unit 
cube. This is achieved by a trivial change of variables 4> — > 4> = <^>/27r an d x — > x, where x 
is defined by 

x (x, v) = x (v) + (1 — x (v)) x, (4.2) 

where x(v) is the lower limit on the x integration. Numerical implementation of the B 
distribution requires all plus distributions be replaced by regular functions, the results of 
which are given in Appendix. |A|. 

The generation of the leading-order, iV-body configuration proceeds as follows: 

1. a leading-order configuration is generated using the standard Herwig++ leading-order 
matrix element generator, providing the Born variables <3?£ with an associated weight 

2. radiative variables are then generated by sampling B ($b), parametrized in terms 
of the 'unit-cube' variables x, v, (j>, using the Auto-Compensating Divide-and-Conquer 
(ACDC) phase space generator [43], which implements a variant of the VEGAS algo- 
rithm [44]; 

3. the leading-order configuration is accepted with a probability proportional to the 



integrand of Eq. 44 evaluated at {p 2 , y, 



9 Apart from the requirement that the final state be colour neutral, this reweighting is independent of 
the details of the vector boson and any decay it undergoes. 
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4.2 Generation of the hardest emission 



The hardest (highest pr) emission is generated from the iV-body configuration according to 
the modified Sudakov form factor, Eq. 2.3. The integrand in the exponent of the Sudakov 
form factor consists of three different contributions, one for each channel ab=qq, qg, gq. 
The integrands are defined as 

tit ^ \ Rab($B,$R) a S C a b 1 ^ S , v , , 

(* fl , $b) = — g($ B ) — = 2 7t x ( x e, ^e) > (4-3) 

where W a fe is equal to 7i a b without the plus prescription. We use the following parametriza- 
tion of the partonic, vector boson and jet momentum in the hadronic centre-of-mass frame: 

P® = T^Vs ( X ®,Q,Q,+ X ®) > p = (m T coshy, p T sm(f), p T cos (f), m T sinh y) , 

Pe = nVs{x e ,0,0, -x e ) , k = (p T coshy k , -p T sin0, -p T cos <f>, p T sinh y k ) , (4.4) 



where pt is the transverse momentum, wit = y P 2 + Pt and is the rapidity of the 
additional parton. Instead of generating the hardest emission in terms of = {x, v, cp} 
we find it more convenient to make a change of variables to <&' R = {pt, yk, related to 
&r according to 

1 4_£.< 1 -„ )(1 -x)», ^^(^-^(i^l). ,4,) 



The modified Sudakov form factor of Eq.l^ contains a ^-function in pt, however by 
choosing to parametrize the radiative phase space in pt the ^-function simply becomes the 
upper limit on the integral. The modified Sudakov form factor for each channel therefore 
has the form 

A R ab (Pt) = exp J d$ R W ab $ B )) , (4.6) 

where PTm&x is the maximum possible transverse momentum. The full Sudakov form factor, 
Aft(py) is given by the product of A a (pt) for the individual channels. The radiative vari- 
ables (pt, yk) are generated according to Eq.|2l| using the veto algorithm. 10 This procedure 
requires simple bounding functions for each channel. Functions of the form, 

9ab(PT) = ^, (4.7) 

Pt 

are used, with suitable values of K a b for each channel together with an overestimate of the 
limits for the rapidity integral, yk min and 2/fc max - The generation procedure then proceeds as 
follows: 



1. Pt is Set to PTmax! 



10 A good description of this technique can be found in Ref. [1]. 
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2. a new (px, Uk) configuration is generated using two random numbers according to 11 

p T =(—-—- -In 7^ , (4.8a) 

\p T K ab (y fcmax - y kmi J J 

Vk = Vk min + K (yfc max - yk m J ; (4.8b) 

3. if pt < PTmin, where PTmin is the minimum allowed pt for the emission, then no 
radiation is generated; 

4. if the generated configuration is outside of the exact phase space boundaries then 
return to step §; 

5. if Wab ($b, 3>h) /dab (pr) > then accept the configuration, otherwise return to step 

!■ 

For this process there are three partonic channels contributing to the radiative corrections, 
this is dealt with by using competition, where a (pt, Uk) configuration is generated, as 
outlined above, for each channel individually and the configuration with the highest px 
accepted. 

We employ a simple prescription [13] to generate the azimuthal angle that allows the 
leptonic correlations to be correctly generated. For the qq channel, the prescription proceeds 
as follows: 

1. momenta are first constructed in the vector boson rest frame; 

2. the Pq, direction is chosen with probability 

(s + i) 2 /({s + i) 2 + (s + u) 2 ), (4.9) 

otherwise the p Q direction is chosen. The momenta are then rotated around the 
chosen direction by a random angle generated uniformly on the interval [0, 2ir] ; 

3. momenta are boosted back to the lab frame such that the rapidity of the vector boson 
is the same as for the iV-body configuration. 

The same procedure is used for the qg and gq initiated channels with the replacements 
s — > t, i — > -ft, u ^ s and s — > u, t — > s, u — ► t, respectively. 

4.3 Truncated and vetoed parton showers 

Before describing how the radiative events are further showered, we need to recall some 
details of the Herwig++ parton shower algorithm. This is described in more detail in 
Refs. [18,39]. The shower starts at a scale given by the colour structure of the hard 
scattering process and evolves down in the evolution variable q by the emission of partons 
in 1 — ► 2 branching processes. Finally, the set of scales, q, momentum fractions, z, and 
azimuthal angles, (j), which describe these branchings, are used to construct the momenta of 

11 1Z refers to a random number in the interval [0, 1], a different random number is generated each time. 
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all the partons radiated in the parton shower. The Herwig++ approach generally requires 
some reshuffling of these momenta after the generation of the parton showers in order to 
ensure global energy-momentum conservation. 



The iV+l-body states generated as described in Sect. [L2J are first interpreted as a single 
standard Herwig++ shower emission, from the iV-body configuration, described by the 
shower branching variables (q~h, Zh, <fih)- The complete POWHEG shower is then performed 
as a single Herwig++ shower modified by certain conditions which allows a simple but 
complete implementation of the truncated shower. 

The shower algorithm proceeds as follows: 

1. the truncated shower evolves from the hard scale, determined by the colour structure 
of the iV-body process, to the hardest emission scale q~h such that the px is less than 
that of the hardest emission, the radiation is angular-ordered and branchings do not 
change the flavour of the emitting parton; 

2. the hardest emission is forced with shower variables (q~h, Zh, <f>h)', 

3. the shower is allowed to evolve down to the hadronization scale with the addition of 
a transverse momentum veto on radiation above px h ■ 

The implementation described above requires that the hardest emission, generated as de- 



scribed in Sect. |4.2| , be interpreted as a Herwig++ shower emission that is forced when we 
evolve to the associated scale in the parton shower. In order to do this, we need to find a 
mapping from the N + 1 momenta, describing the hardest emission, to the shower variables 
(<27i, Zh, 4>h) and an iV-body configuration. This equates to undoing the momentum recon- 
struction procedure used in the Herwig++ shower. The reconstruction procedure consists 
of two steps [39]. First, the momenta of partons at each step of the shower are constructed, 
in the centre-of-mass frame of the hadronic collision, recursively from the shower variables. 
Second, boosts are applied to each jet individually which ensures global momentum con- 
servation. The reconstruction process is different for initial- and final-state radiation, here 
we will only consider the initial-state case which is relevant for the Drell Yan process. 
In Herwig++ the momenta of the partons, g,, in a jet are given by 

qi = aip + fan + qj_i, (4.10) 

where for initial-state radiation the reference vectors p and n are given by the hadronic 
momenta of the beam particles and Pq and q±i is the transverse momentum with respect 
to the beam axis. 

For initial-state radiation we use a backward evolution algorithm which starts from the 
hard process and evolves to lower evolution scales backwards towards the incoming hadron 
by the emission of time-like partons. The reconstruction of the initial-state jet starts from 
the last initial-state parton produced by the backward evolution algorithm with momentum 
calculated from the fraction of the beam momentum it carries. The momentum of the time- 
like daughter of this parton is reconstructed as described in Ref. [39]. The momentum of 
the space-like daughter is then given by momentum conservation. This process is iterated 
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for each initial-state branching eventually giving the momentum of the space-like progenitor 
parton, colliding in the hard process. 

The momenta of the two progenitor partons are then reshuffled such that mass and 
rapidity of the partonic centre-of-mass system is conserved. Under this reshuffling the 
progenitor momenta are transformed according to 



q®^q'® = a® k ®p® + t^p® + ff±@- (4.ii ■ 

K ® 



The reshuffling parameters, fc® and fc e , are found by solving equations requiring conserva- 
tion of mass and rapidity and hence the associated Lorentz transform is obtained. 



The basis vectors are the hadronic beam momenta and the a@ parameters in Eq. |4.11 



are simply the Born partonic momentum fractions, given in Eq, |3.2| . The reshuffling pa- 
rameters, can therefore be calculated from the momenta of the shuffled progenitors, 
c/o@- Decomposing these momenta into their Sudakov parameters the momentum shuffling 
parameters are simply 

k 9 = ^, k e = ^e, (4.12) 

Xq Xq 

where a' 0@ refers to the a parameters in the Sudakov decomposition of the shuffled progeni- 
tors. The inverse of the Lorentz boosts implementing the reshuffling can then be calculated 
and applied to each momentum, yielding the unshuffled momenta qi. These are then decom- 
posed yielding their Sudakov parameters, the shower variables parametrizing the branching 
can then be determined. The momentum fraction is given by 

z=^, (4.13) 

where at is the Sudakov parameter for the space-like parton entering the hard process 
and a-j the Sudakov parameter of the initial-state parent parton. In this simple case the 
transverse momentum is simply equal to that of the off-shell space-like parton initiating the 
leading-order hard process, or equivalently, its outgoing, time-like, sister parton. The scale 
of the branching is defined in terms of the pr and light-cone momentum fraction z, as 

f = Z -^±f, (4.14) 

where Q g is the constituent gluon mass, the infrared regulator of the Herwig++ parton 
shower. 

Using this approach we can calculate the shower variables (q~h, Zh-,4> z ) for the hardest 
emission which allows us to generate the truncated and vetoed showers. 

5. Results 

As a check of the calculation of the next-to-leading order differential cross section, distri- 
butions of the vector boson rapidity produced by the POWHEG implementation and the 
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a) Z Tevatron b) W + Tevatron 




Figure 1: Comparisons of da/dy for the POWHEG implementation and MCFM [45] for Z and W + 
production at the Tevatron (y/s = 2TeV) and the LHC (y/s = 14TeV). 



NLO program MCFM [45] were compared. Fig.[| shows distributions for 7/Z and W + pro- 
duction at the Tevatron (proton-antiproton at yfs = 2 TeV) and the LHC (proton-proton at 
y/s = 14 TeV). In all cases the total cross sections from MCFM and our POWHEG implemen- 
tation agreed to within 0.5 %. The distribution of the rapidity of the lepton produced in the 
j/Z and W decay is shown in Fig. |2] and is also in good agreement. For both Herwig++ and 
MCFM in this comparison, the parton density functions used were the MRST2001 NLO [46] 
set with the LHAPDF interface [47]. 

In Figs.||-||, distributions from the Drell-Yan POWHEG implementations for the ra- 
pidity and transverse momentum of the vector bosons are compared to Tevatron data. The 
middle and bottom panels in each of these plots shows the (Theory — Data) /Data and x 
values for each bin. In Fig. || the rapidity distribution of j/Z bosons of mass 71-111 GeV 
is compared to DO Run II data [48]. Fig. |I] shows the transverse momentum distribu- 
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Figure 2: The rapidity of a) the electron in Z and b) the positron in W + production at the 
Tevatron including the leptonic decay of the gauge boson for the POWHEG implementation and 
MCFM [45] at the Tevatron (y/s = 2TeV). 

tion of j/Z bosons of mass 66-116 GeV compared to CDF Run I data [49]. Fig. || shows 
the transverse momentum distribution of ~f/Z bosons of mass 40-200 GeV compared to DO 
Run II data [50]. Fig. || shows the transverse momentum distribution of W bosons com- 
pared to Run I DO data [51]. In addition to the results from our implementation of the 
POWHEG method the results from Herwig++ including a matrix element correction and 
MCQNLO [19-24] are shown. The predicted W and Z px distributions at the LHC are 
shown in Fig.|7[ 

The Herwig++ results were generated using an intrinsic px of 2.2 GeV which was ob- 
tained by fitting to the Run I W and Z px distributions [39]. The POWHEG results used 
the same intrinsic pt as for Herwig++ and a minimum px of 2 GeV for the hardest emission. 
The MCONLO and HERWIG results were generated using an intrinsic px of 1.6 GeV from a 
fit to DO data [52]. 

The leading-order parton distribution functions of [46] were used for the Herwig-I— (- 
result and the central value of the NLO parton distributions of [53] for the POWHEG and 
MCONLO results. 

All the approaches give good agreement for the rapidity of the Z boson however they 
differ in the description of the px spectrum of the gauge boson. The chi squared per 
degree of freedom for the various px spectra and approaches are given in Table ffl. All the 
approaches are in good agreement with the Run I data from CDF and DO for the px of the 
Z and W. However, with the exception of the results of the FORTRAN HERWIG program 
including a matrix element correction, which gave the worst agreement with the Run I Z 
data, all the results are below the new DO Z px data at high transverse momentum. 
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Figure 3: Rapidity distribution for Z production compared to DO Run II Tevatron data [48]. The 
solid line shows the prediction of our POWHEG implementation, the dotted line is the prediction 
of MC@NLO and the dashed line is the default Herwig++ result. 

There is a common trend, at both the Tevatron and LHC energies, that the matrix 
element correction gives the largest result at large pt, followed by the POWHEG approach 
with MC@NLO giving the lowest value. This is due to the treatment of the hardest emission 
in the different approaches. In the MCONLO method the result at large px is the leading- 
order matrix element for the production of a vector boson and a hard QCD jet. However in 
this region, as we are normalising to the total cross section, the matrix element correction 
result is essentially the matrix element for vector boson plus jet production multiplied by 
the K-factor 12 giving a larger result. In the large px region the POWHEG result, because 
the real-emission matrix element is exponentiated, is the real-emission matrix element mul- 
tiplied by the B function, which results in a K-factor-like correction, and the Sudakov form 
factor which causes the result to be slightly smaller than the default Herwig++ result. The 
POWHEG result has the significant advantage that rather than using a global rescaling 
of the cross section to get the NLO normalization, which can lead to a poor description 
of observables, such as the boson rapidity, which are non-zero at leading order the NLO 
correction is calculated for each momentum configuration. 

In an ideal world we would like to use the NLO result for vector boson production 
in association with a hard jet to describe the high pt region, however incorporating this 
result into a Monte Carlo simulation is not currently feasible, and therefore the POWHEG 

12 The K-factor here is the ratio of the NLO cross section for inclusive vector boson production divided 
by the leading-order cross section. 
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Figure 4: Transverse momentum distribution for Z production compared to CDF Run I Tevatron 
data [49]. The solid line shows the prediction of our POWHEG implementation, the dotted line is 
the prediction of MC@NL0 and the dashed line is the default Herwig++ result. The inset shows an 
expanded view of the low pr region. 



Approach 


Data Set 




DO W pt 


CDF Z p T 


DO Z p T 




All 


p T > 30 GeV 


All 


p T > 30 GeV 


All 


p T > 30 GeV 


MC@NL0 


0.51 


0.82 


0.70 


0.96 


7.2 


13.9 


Herwig++ 


0.67 


0.42 


0.89 


0.61 


5.1 


7.0 


POWHEG 


0.54 


0.33 


1.99 


1.00 


5.3 


6.9 


HERWIG 


0.69 


1.08 


2.45 


4.47 


2.0 


1.9 



Table 1: Chi squared per degree of freedom for MC@NLO, Herwig++, our implementation of the 
POWHEG method in Herwig++ and FORTRAN HERWIG compared to Tevatron vector boson pr 
data. The chi-squared values are calculated for the shapes of the distributions, i.e. normalizing 
them to unity. In order to compare the high pr region and minimise the effect of tuning the intrinsic 
transverse momentum the chi squared per degree of freedom is given for both the full pt region 
and only for the data points with pr > 30 GeV. 

or matrix element correction methods which basically assume the correction for inclusive 
production is the same as for vector boson production in association with a jet at least have 
the advantage of including a correction which improves agreement with data. 

In general all the results lie below the DO Run II Z pt data between 50 and 100 GeV 
which results in the relatively poor chi squared, however in general the POWHEG approach 
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Figure 5: Transverse momentum distribution for Z production compared to DO Run II Tevatron 
data [50]. The solid line shows the prediction of our POWHEG implementation, the dotted line is 
the prediction of MC@NL0 and the dashed line is the default Herwig++ result. The inset shows an 
expanded view of the low pr region. 

gives comparable results to the other state-of-the-art techniques. The effect of varying the 
scale used for the parton distributions and as between 0.5s and 2s for the B term and 
between 0.5(Mg + p T ) and 2(M| + p T ) for the hardest emission is shown in Fig.||. While 
this variation moves the POWHEG result close to the data it still is below the experimental 
result in the intermediate px region. 

The effect of the truncated shower is illustrated in Fig.|| which shows the small pr 
region of the transverse momentum distribution for W and Z production compared to DO 
and CDF data. In this region where the highest px emission is at a small scale and there 
is often a large region for the evolution of the truncated shower it has the largest effect. 
However the effect is relatively small at least for the transverse momentum distribution, 
equivalent to a small change in the intrinsic transverse momentum. We would expect to 
see a larger effect in the distributions of jets in the event, in particular the second hardest 
jet, which we will study in the future. 

6. Conclusion 

The POWHEG NLO matching prescription has been implemented in the Herwig++ Monte 
Carlo event generator for Drell-Yan vector boson production. A full treatment of the 
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Figure 6: Transverse Momentum distribution for W production compared to DO Run I data [51]. 
The solid line shows the prediction of our POWHEG implementation, the dotted line is the prediction 
of MC@NL0 and the dashed line is the default Herwig++ result. The inset shows an expanded view 
of the low pt region. 

truncated shower, which is required to produce wide angle, soft radiation in angular- ordered 
parton showers, is included for the first time. 

The implementation gives a good description of data, on a similar level to the matrix 
element correction methods and better than MCONLO. It will be available in a forthcoming 
release of Herwig++. 

The technique we have used to implement the POWHEG approach, by interpreting 
the hard emission in terms of the variables used to generate the parton shower is very 
powerful, and has many other applications in matching approaches, for example for the 
CKKW approach [32,33], which we will explore in the future. 
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a) W Production b) Z Production 
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Figure 7: The px distributions of a) W and b) Z bosons at the LHC. The solid line shows the 
prediction of our POWHEG implementation, the dotted line is the prediction of MC@NL0 and the 
dashed line is the default Herwig++ result. 
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Figure 8: Transverse Momentum distribution for Z production compared to DO Run II data [50]. 
The band shows the effect of varying the scale used for the parton distributions and as between 
0.5s and 2s for the B term and between 0.5(M§ + p T ) and 2(M§ + p T ) for the hardest emission. 
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Figure 9: Transverse Momentum distribution for a) W production compared to DO Run I data [51] 
and b) Z production compared to CDF Run I Tevatron data [49]. The solid line includes the 
truncated shower whereas the dashed line does not. 

Note added in Proof 

While we were in the final stages of completing this paper another paper on the same topic 
was submitted to the arXiv [54]. 



A. Plus Distributions 



In order to implement the collinear (C a b) terms in the real-emission contributions to B ($b) 
the following relations are required 



IL) dx 



(v) aX (1-x) 



and 



J^- + = j, 1 dx (i - x ( V) ) + Wr ln ( x - * 



(A.l) 



fi dx (1 - x (v)) [(/ (x\x, v))-f(x (1, «))) ( 'TS^ ) + ln2 C 1 " * ^ 

(A.2) 

with x defined in Eq. 4. 10| and v £ [0,1]. For the hard (H a b) contribution to the real 
radiation components in B ($b) 

Jo dv dx f ( x > v ) (T^)7 ( (T^)7 + i) 

- ^rh, f 1 ^? 1 f /feMMM!M !) I /(a(g,«),«)-/(l,«)-/(g(g,0),0)+/(l,0) 

~~ Jo ut7 Jo ux i-a ^ i-o u 



+ Jo du Jo dx 



y 1—V ' V 

l,^)ln(l-5i( t) ))-/(l,l)ln(l-5:(l)) f(l,v) ln(l-x(v))-f(l,0) ln(l-z(0)) \ 
1—v v J 



we 



where in the last line of Eq. A. 3 , 

are derived, in different variables, in Ref. [29] 



(A.3) 



have introduced the identity as f Q dx. Similar relations 
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